#### This code builds the doble branched covers of all alternating links in the HT link table.
#### If they homology of the DBC D over the link L is cyclic, the name of L, the order of the homology of D and the crossing number of L are saved in a csv-file.
#### On my computer it takes about 45 min.

import time
import pandas
import snappy

def double_branched_cover(link):
    """
    Returns the double branched cover of the link.
    """
    L=link.copy()
    for i in range(L.num_cusps()):
        L.dehn_fill((2,0),i)
    for cov in L.covers(2):
        if (2.0, 0.0) not in cov.cusp_info('filling'):
            return cov

start_time = time.time()
pandas.DataFrame({'knot': [], 'homology': [], 'crossings': []}).to_csv('DBChomologies.csv')

for L in HTLinkExteriors(alternating=True):
    D=double_branched_cover(L)
    if len(D.homology())==1:
        df=pandas.DataFrame({'knot':L.name(), 'homology':D.homology().coefficients, 'crossings':len(L.link().crossings)})
        df.to_csv('DBChomologies.csv', mode='a', header=False)
print("--- Time taken: %s seconds ---" % (time.time() - start_time))

# Remark: If one is just interested in the order of the homologies one can also use SnapPy inside sage and compute the determinant of the links which is the order of the homology of the DBC of that link.

### It is also possible to directly save the triangulations of the DBCs. For that use M._to_string() to get a string that describes the manifold. To load the manifold again use M(string).

########################################################################################################
######### We also do the same for the non alternating links. This is needed for the search of dbc slopes.

start_time = time.time()
pandas.DataFrame({'knot': [], 'homology': [], 'crossings': []}).to_csv('DBChomologies_non_alternating.csv')

for L in HTLinkExteriors(alternating=False):
    D=double_branched_cover(L)
    if len(D.homology())==1:
        df=pandas.DataFrame({'knot':L.name(), 'homology':D.homology().coefficients, 'crossings':len(L.link().crossings)})
        df.to_csv('DBChomologies_non_alternating.csv', mode='a', header=False)
print("--- Time taken: %s seconds ---" % (time.time() - start_time))


### Both searchs take around 1 hour.

#####################################################################################
### Next we build data bases for surgery diagrams as well.

def gcd(a, b):
    while b != 0:
        a, b = b, a % b
    return a

def coprime(a, b):
    return gcd(a, b) == 1

slope_set=[]
for q in range (1,5):
    for p in range (2,11):
        if coprime(p,q)==True:
            slope_set.append((p,q))
            slope_set.append((-p,q))
            
            
def double_branched_cover(link):
    """
    Returns the double branched covers of the link. This works also for links in a more general manifold. 
    Note that a knot in a more general manifold may have more than one double branched cover 
    (or no double branched cover at all if the knot represents a primitive element in homology).
    This function will return the complete list of all double branched covers of the link.
    """
    L=link.copy()
    for i in range(L.num_cusps()):
        if L.cusp_info(i).filling==(0.0, 0.0):
            L.dehn_fill((2,0),i)
    return [cov for cov in L.covers(2) if (2.0, 0.0) not in cov.cusp_info('filling')]
    
    
    


start_time = time.time()
pandas.DataFrame({'knot': [], 'cusp': [], 'filling': [], 'homology': [], 'crossings': []}).to_csv('DBChomologies_one_filling.csv')

for i in range(1,13):
    for L in snappy.HTLinkExteriors(crossings=i,knots_vs_links='links'):
        for i in range(0,L.num_cusps()):
            for r in slope_set: 
                L.dehn_fill(r,i)
                DBC=double_branched_cover(L)
                for D in DBC:
                   if len(D.homology())==1:
                       df=pandas.DataFrame({'knot':L.name(), 'cusp':[i], 'filling':[r], 'homology':D.homology().coefficients, 'crossings':len(L.link().crossings)})
                       df.to_csv('DBChomologies_one_filling.csv', mode='a', header=False) 
                L.dehn_fill((0,0),i)
print("--- Time taken: %s seconds ---" % (time.time() - start_time))  


#### This takes around 2 hours. The same for more crossings and more fillings (which takes another 2,5 days). In total we have constructed around 7.7 million surgery diagrams.

start_time = time.time()
for i in range(13,15):
    for L in snappy.HTLinkExteriors(crossings=i,knots_vs_links='links'):
        for i in range(0,L.num_cusps()):
            for r in slope_set: 
                L.dehn_fill(r,i)
                DBC=double_branched_cover(L)
                for D in DBC:
                   if len(D.homology())==1:
                       df=pandas.DataFrame({'knot':L.name(), 'cusp':[i], 'filling':[r], 'homology':D.homology().coefficients, 'crossings':len(L.link().crossings)})
                       df.to_csv('DBChomologies_one_filling.csv', mode='a', header=False) 
                L.dehn_fill((0,0),i)
print("--- Time taken: %s seconds ---" % (time.time() - start_time)) 

###  and with more slopes:

slope_set=[]
for q in range (1,3):
    for p in range (1,5):
        if coprime(p,q)==True:
            slope_set.append((p,q))
            slope_set.append((-p,q))

start_time = time.time()
pandas.DataFrame({'knot': [], 'cusp1': [], 'filling1': [], 'cusp2': [], 'filling2': [], 'homology': [], 'crossings': []}).to_csv('DBChomologies_two_filling.csv')

for i in range(1,15):
    for L in snappy.HTLinkExteriors(crossings=i,knots_vs_links='links'):
        if L.num_cusps()>2:
            for i in range(0,L.num_cusps()):
                for r in slope_set: 
                    L.dehn_fill(r,i)
                    F=L.filled_triangulation([i])
                    for j in range(0,F.num_cusps()):
                        for s in slope_set:
                            F.dehn_fill(s,j)
                            DBC=double_branched_cover(F)
                            for D in DBC:
                                if len(D.homology())==1:
                                    df=pandas.DataFrame({'knot':L.name(), 'cusp1':[i], 'filling1':[r], 'cusp2':[j], 'filling2':[s], 'homology':D.homology().coefficients, 'crossings':len(L.link().crossings)})
                                    df.to_csv('DBChomologies_one_filling.csv', mode='a', header=False) 
                L.dehn_fill((0,0),i)
print("--- Time taken: %s seconds ---" % (time.time() - start_time))  
              
